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Abstract 

In this paper we discuss techniques, which lead to a significant improvement of the efficiency 
of the Monte Carlo integration, when one-loop QCD amplitudes are calculated numerically 
with the help of the subtraction method and contour deformation. The techniques discussed 
are: holomorphic and non-holomorphic division into sub-channels, optimisation of the inte- 
gration contour, improvement of the ultraviolet subtraction terms, importance sampling and 
antithetic variates in loop momentum space, recurrence relations. 



1 Introduction 



In a recent letter [HI we reported on the next-to-leading order (NLO) results in the leading-colour 
approximation for five, six and seven jets in electron-positron annihilation. This was the first time 
a physical observable depending on a one-loop eight-point function has been calculated. The 
calculation has been performed with a method, which uses numerical Monte Carlo integration 
for the computation of the one-loop amplitude. To render this part finite, the subtraction method 
El-El and a contour deformation [|2ll3ll6T-fTT| have been used. Also, it has been the first time 
ever that the numerical method has been applied to a cutting-edge process. Competing groups 
use either unitarity techniques [fl2~l - l2"3l or Feynman diagram techniques ll24ll . For the unitarity 
method or the Feynman diagram method a variety of packages are available lf25l - [33Tl . 

Within all approaches for multi-parton NLO calculations, the computation of the virtual cor- 
rections is the most challenging part. The basic principles how the virtual corrections are calcu- 
lated within the numerical method have already been discussed in the literature El-flTI. What is 
missing in the literature is information on how this can be done efficiently. Efficiency is a cru- 
cial ingredient to apply the method to high multiplicity processes involving seven or eight point 
functions. 

In this paper we provide detailed information how the numerical method for the computation 
of the virtual corrections can be implemented in an efficient way. There are two important parts: 
The first essential part is the reduction of the statistical error of the Monte Carlo integration for 
the virtual part. We employ several techniques to achieve this goal. The techniques used are: 
holomorphic and non-holomorphic division into sub-channels, optimisation of the integration 
contour, improvement of the ultraviolet subtraction terms as well as importance sampling and 
antithetic variates in loop momentum space. The numerical method uses contour deformation 
and subtraction terms to avoid singularities. In the vicinity of a singularity we use whenever 
possible contour deformation to by-pass the singularity. In the case where this is not possible, 
because the contour is pinched, and if the singularity is not integrable then there is a subtraction 
term for it. This occurs in the soft and collinear regions. We expect that the regions close to soft 
and collinear singularities are delicate with respect to the Monte Carlo integration and will give 
a significant contribution to the overall Monte Carlo integration error. This is indeed the case 
and a solution of this problem belongs in the category of solutions for the "known unknowns". 
It turned out that we also needed a solution from the category of the "unknown unknowns". In a 
first attempt we started with an integrand which falls off like \k\~ 5 for \k\ — > °°. Formally this is 
sufficient to ensure ultraviolet finiteness. However, it turned out that the ultraviolet region gave 
a large contribution to the overall Monte Carlo error. A significant part of this paper is devoted 
on the improvement of the ultraviolet behaviour. 

The second part is the efficient calculation of the integrand. To achieve this goal, recurrence 
relations are used. We show how the integrand of the bare one-loop amplitude and the subtraction 
terms can be calculated efficiently. In essence the calculation of the integrand of the bare one- 
loop amplitude reduces to a tree-like recurrence relation, once the loop has been cut open in a 
single place. 

This paper is organised as follows: In the next section we briefly review the basic principles 
of the numerical method for the computation of the virtual corrections. In section [3] we discuss 
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in detail the techniques which we use to reduce the statistical Monte Carlo error. Section @] is 
devoted to recurrence relations. Finally, section |5]contains our conclusions. 



2 Numerical calculation of one-loop amplitudes 

In this section we briefly outline the method for the numerical computation of NLO corrections. 
More details can be found in ||2). For concreteness we discuss the case of electron-positron 
annihilation. However all methods can equally well be applied to hadron-hadron collisions and 
deep-inelastic scattering. 



2.1 Overview of the method 

In electron-positron annihilation the contributions at leading and next-to-leading order for an 
infrared-safe observable O are given as 

LO f ^ j — -B /^>\NLO n , / r\ j^V 



(<2} LU = J O n do a , (<9) 1NLU = J O n+l do K + j O n da\ (1) 

n n+l n 

Here a rather condensed notation is used. do B denotes the Born contribution, whose matrix 
elements are given by the square of the Born amplitudes with (n + 2) particles. Similar, do R 
denotes the real emission contribution, whose matrix elements are given by the square of the Born 
amplitudes with (n + 3) particles. Jo v gives the virtual contribution, whose matrix elements 
are given by the interference term of the one-loop amplitude with (n + 2) particles, with the 
corresponding Born amplitude. Taken separately, the individual contributions at next-to-leading 
order are divergent and only their sum is finite. In order to render the real emission contribution 
finite, such that the phase space integration can be performed by Monte Carlo methods, one adds 
and subtracts a suitably chosen piece [|34| - |371 : 



(0> NLO = J (O n+l da R - O n do A ) + J I O n da y + O n J do A j 

n+l n \ I ) 



(2) 



The term {O n +\ do R — O n do A ) in the first bracket is by construction integrable over the (n + 1)- 
particle phase space and can be evaluated numerically. The subtraction term can be integrated 
analytically over the unresolved one-particle phase space. In a compact notation the result of this 
integration is often written as 



/ 



do A = I®da B . (3) 
l 

The notation <g> indicates that colour correlations due to the colour charge operators still remain. 
The virtual contribution do y is given by 

da y = 2Re ^ (0) *J? (1) )^, 2 . (4) 
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jftW denotes the renormalised one-loop amplitude. It is related to the bare amplitude by 



1) 



-^bare + ^CT ■ 



(5) 



,(1) 



denotes the ultraviolet counterterm from renormalisation. The bare one-loop amplitude 
involves the loop integration 



d D k (i) 

D ^bare ' 



(27!) 



(6) 



where denotes the integrand of the bare one-loop amplitude. Within our approach we 
extend the subtraction method to the integration over the virtual particles circulating in the loop. 
To this aim we rewrite eq. © as 

_ (JX) _aW _ a^A + /V 1 ) + a« + $?« + m 

-Tjare^-^CT _ ^-^bare ^soft ^coll "tJV y ^ ^"r^T ^ soft ' coll ' "tJV J " 

The subtraction terms ^LqI, and J^y-y are chosen such that they match locally the singular 
behaviour of the integrand of in D dimensions. The first bracket in eq. © can therefore 
be integrated numerically in four dimensions. The term J^L approximates the soft singulari- 
ties, J? c ^jj approximates the collinear singularities and the term J3^y approximates the ultraviolet 
singularities. These subtraction terms have a local form similar to eq. ©: 



soft 



d D k (i) 

\D ^soft' 
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(1) 



r (l) 

\£> ycoii' 



/ 



r (l) 

d yuv 



(271)° r7sotp "" co11 7 (27r) Dr7co11 ' " aJV 7 (2tc) 
The contribution from the terms in the first bracket of eq. (J7} can be written as 



(8) 



/ 



2 Re 



l-^bare -^soft 



coll 



2 Re 



G W ■ 

^bare 



r (l) 
Vsoft 



ycoii 



r (i) 
yuv 



On +0(8). 



(9) 



The integral on the right-hand side is finite. Note that the integration over the loop momentum k 
can be done together with the phase space integration in a single Monte Carlo integration. The 
building blocks of the subtraction terms are process-independent. When adding them back, we 
integrate analytically over the loop momentum k. The result can be written as 



2 Re 



-^CT + Aoft + -^coll + AJV 



L<g>dcr 



(10) 



The insertion operator L contains the explicit poles in the dimensional regularisation parameter 
related to the infrared singularities of the one-loop amplitude. These poles cancel when combined 
with the insertion operator I: 



(I + L) 0do t 



finite. 



(11) 
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For massless particles we have 



I + L 



2ft 



Re 



EET f T; 

1 j& 



In 



7T- 



6(2^0 



2 fj z 



+ 0(8). 



(12) 



In this formula we denote by /j the renormalisation scale and by puy the scale used in the ultra- 
violet subtraction terms. Of course it is possible to set /j = fj\jw, but it will be advantageous to 
keep this two scales different. In fact, we will choose /j^y purely imaginary. The colour charge 
operator for particle i is denoted by T,. We further have 



Cp, 



Kq — Kq 



ft 

~6 



Cp, 



K„ 



67 _T? 
18 ~~6~ 



C A 



10. 



T R N f . 



(13) 



Po is the first coefficient of the QCD |3-function and given by 

Po=—C A --T R N f . 



The colour factors are as usual 



Nl-\ 

C A =N C , C F = —-.. 



Tr 



1 

2' 



(14) 



(15) 



We therefore have to evaluate three contributions to get the next-to-leading order correction: The 
real emission contribution 



(^)reai° 



J (O n+ ido R -O n da A ) 



(16) 



B+l 



the insertion term 



/<9\ NL °. 

\^/ insertion 



J 0„(I + L)®da E 



(17) 



and the virtual contribution 

<0>S£d = z/^Re 



I %are ^soft 



ycoii 



uv 



O n . (18) 



All three contributions are by construction finite. In this paper we focus on the virtual contribu- 



5 



2.2 Colour decomposition and kinematics 

It is convenient to decompose a full one-loop QCD amplitude into primitive amplitudes: 

A® = Z,CjAf\ (19) 

j 

The colour structures are denoted by Cj while the primitive amplitudes are denoted by Aj. 
In the colour-flow basis [|38l - l40l the colour structures are linear combinations of monomials in 
Kronecker 8,-/s. A primitive amplitudes is defined as a colour- stripped gauge-invariant set of 
Feynman diagrams with a fixed cyclic ordering of the external partons and a definite routing of 
the external fermion lines through the diagram [@T|- 

It is simpler to work with primitive one-loop amplitudes instead of a full one-loop amplitude. 
Our method exploits the fact that primitive one-loop amplitudes have a fixed cyclic ordering 
of the external legs and that they are gauge-invariant. The first point ensures that there are 
at maximum n different loop propagators in the problem, where n is the number of external 
legs, while the second property of gauge invariance is crucial for the proof of the method. We 
therefore consider in the following just a single primitive one-loop amplitude, which we denote 
by A^ l \ while keeping in mind that the full one-loop amplitude is just the sum of several primitive 
amplitudes multiplied by colour structures. The full one-loop amplitude is obtained by summing 
over all primitive amplitudes. 

We introduce some notation related to the primitive amplitude A^\ Since the cyclic ordering 
of the external particles is fixed, there are only n different propagators occurring in the loop 
integral. We label the external momenta clockwise by p\, p2, p„ and define qi = p\ + P2 + 
... + pi. The loop momenta are 

j 

1=1 

For convenience we set 

ko = k„ and qo = q n - (21) 

Due to momentum conservation we actually have 

qo = q n = 0. (22) 

Nevertheless we will use go (or this makes the formulae more symmetric with respect to 
the indices. In electron-positron annihilation we can take pi, p2, p n -2 t0 be the final state 
momenta and p n -\ and p n to be the (negative) of the initial state momenta. The two leptons cou- 
ple only through a photon or through a Z-boson to the quark line. This implies that in primitive 
amplitudes related to electron-positron annihilation only {n — 1) loop propagators are present. 
Phrased differently, in electron-positron annihilation the propagator corresponding to the loop 
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momentum k n -\ is absent. We can write the bare primitive one-loop amplitude in Feynman 
gauge as 



A ^ e ~ / (2%) D ^ G baL " ^are (*) \j fe2 _ ^2 + f g " ( 23 ) 

G^2e is the integrand of the bare one-loop amplitude. Pbare(^) is a polynomial in the loop mo- 
mentum k. The +z'8-prescription instructs us to deform - if possible - the integration contour into 
the complex plane to avoid the poles at kf = mf. In the following we do not write the +z'8- term 
explicitly. 



2.3 The subtraction terms 



At the level of primitive amplitudes we denote the subtraction terms for the integrand, corre- 
sponding to the ones appearing in eq. (fT8l ) by 



(i) 

'soft' 



(1) 

'coll' 



'UV 



(24) 



For massless QCD the soft and collinear subtraction terms are given by 



'soft 



'coll 



Pj-Pj+i _ A (o) 



i-, u2 k2k2 .1 

ja g K j-\ K j K j+\ 



Sjguv{lCj-i,tf) S j+ ig V y(k 2 ,k 2 +l ] 



k 2 k 2 



+ ■ 



k 2 k 2 



A (0) 



(25) 



where the sum over j G I g is over all gluon propagators j inside the loop. If we take the subset 
of diagrams which have the gluon loop propagator j and if we remove from each diagram of this 
subset the loop propagator j we obtain a set of tree diagrams. After removing multiple copies of 

identical diagrams this set forms a Born partial amplitude which we denote byA^\ Furthermore, 
Sj = 1 if the external line j corresponds to a quark and Sj = 1/2 if it corresponds to a gluon. 
The function guv ensures that the integration over the loop momentum is ultraviolet finite. The 
function guv must have the properties 

lim guv (kj_ u kj) = 1, lim g UV {kj-utf) = o(j). (26) 

kj-i\\kj J JJ k-+o° J J/ \kj 



A possible choice is [0 



#UV 



k 2 k 2 



[(k-Q) 2 -^} 



2 ]2- 



(27) 



Q is an arbitrary four- vector independent of the loop momentum k and is an arbitrary scale. 
Since these two quantities are arbitrary, there are no restrictions on them, they even may have 
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complex values. We will later choose ,u uv purely imaginary with Im /j uv < 0- This will ensure 
that the denominator of eq. (|27T > does not introduce additional singularities for the contour in- 
tegration. There are many possible choices for the function guv compatible with eq. (l26l) . and 
we will choose an improved form which minimises the Monte Carlo integration error. This will 
be discussed in detail in section |3] The soft and collinear subtraction terms for QCD amplitudes 
with massive partons are also known [|2]-|4). The soft and collinear subtraction terms are formu- 
lated at the amplitude level and are proportional to the corresponding Born amplitudes. Upon 
integration they yield simple analytic results: 



r i2 £ f d D k (1) 1 y 2 ( -2pjp j+ i y* {0) 



with S e = (47t) e exp(— £Je) the typical volume factor in dimensional regularisation, where Je 
denotes the Euler-Mascheroni constant, /j denotes the renormalisation scale in dimensional reg- 
ularisation and £ is defined through D = 4 — 2e. 

The ultraviolet subtraction terms correspond to propagator and vertex corrections. The sub- 
traction terms are obtained by expanding the relevant loop propagators around a new ultraviolet 
propagator (k — [J%y)~ , where k = k— Q: For a single propagator we have 

1 1 2k-(p-Q) {p-Q) 2 -m 2 +nl w [2k-(p-Q)] 2 

' "i - 



+ o(4^)- (29) 



1*1 



We can always add finite terms to the subtraction terms. For the ultraviolet subtraction terms we 
choose the finite terms such that the finite parts of the integrated ultraviolet subtraction terms are 
independent of Q and proportional to the pole part, with the same constant of proportionality. 
The integrated ultraviolet subtraction terms have the form 

1- In 40 +0(8), (30) 
£ n 2 J 

where c depends on the type of the subtraction term. This ensures that the sum of all integrated 
UV subtraction terms is again proportional to a tree-level amplitude. 
Putting everything together we have 

r (i) _ r (i) _ r (i) _ r (i) _ R ( k ) 

Uk„™ *~J„ n f, U„„ii <Jt 



bare soft coll UV ;3 _2 

n 
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P(k) and P\jv(k) are, polynomials in the loop momentum k. G h J re and G so f t defined as in eq. (1251) 
contribute only to the first term of the r.h.s of eq. (|3TT i. while G^y contributes only to the second 
term on the r.h.s of eq. (13TT ). On the other hand, G^J n contributes to both terms on the r.h.s of 
eq. CD). 



2.4 The contour deformation 



Having a complete list of ultraviolet and infrared subtraction terms at hand, we can ensure that 
the integration over the loop momentum gives a finite result and can therefore be performed 
in four dimensions. However, this does not yet imply that we can safely integrate each of the 
four components of the loop momentum kf 1 from minus infinity to plus infinity along the real 
axis. There is still the possibility that some of the loop propagators go on-shell for real values 
of the loop momentum. If the contour is not pinched this is harmless, as we may escape into 
the complex plane in a direction indicated by Feynman's +/8-prescription. However, it implies 
that the integration should be done over a region of real dimension 4 in the complex space C 4 . 
If the contour is pinched then the singularity is integrable when the integration is done over the 
loop momentum space and the phase space. This is the case because either the singularity in 
the bare one-loop amplitude is integrable by itself, or - if not - there is a subtraction term for 
it. Let us look again at eq. fUTh . Choosing /j^y sufficiently large on the negative imaginary axis 
ensures that the integration contour stays always away from the poles defined by k 2 — /u^y = 0. 
Therefore we have to consider for the contour deformation only the poles defined by k 2 —m 2 = 0. 
The integral which we will have to consider is given by 



d 4 k 



r (i) _ r m 

^bare ^soft 



'coll 



r UV 



d 4 k 



R(k) 



(271)4-2 
.7=0 



m' 



(32) 



where R(k) is a rational function of the loop momentum which has only poles at k — /Jyv = 
and the integration is over a complex contour in order to avoid whenever possible the poles of 
the propagators shown explicitly in eq. (l32l . 

Let us first look for conditions how the integration contour has to be chosen. We set 



k = k + iK(k) , 
where is real. After this deformation our integral equals 



/ 



dVc_ 



dk? 



R(k(k)) 



n-2 



(33) 



(34) 



LI (k 2 -m 2 -K 2 + 2ikj-K 
j=o \ 



To match Feynman's +?'8-prescription we have to construct the deformation vector K such that 

kj -k>0, (35) 



kl-m 2 
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and the equal sign applies only if the contour is pinched. Integration contours, which fulfil these 
requirements are called allowed contours. Among all allowed contours we would like to pick 
one, 

• which we can construct algorithmically in a process-independent way and 

• which leads to a small Monte Carlo integration error. 

The first requirement is clearly needed for a process-independent method, the second requirement 
is essential for the efficiency of the method. We have implemented and tested two algorithms 
for the contour deformation. The first one uses an auxiliary Feynman parametrisation [|2l[8), 
where also the Feynman parameters are deformed into the complex Feynman parameter space. 
The Monte Carlo integration is then over the phase space of the final state particles, the loop 
momentum space and the Feynman parameter space. The second algorithm follows ref. [0 and 
works directly in loop momentum space. In this case the Monte Carlo integration is only over 
the final state phase space and the loop momentum space. 

The advantages and disadvantages of the two algorithms are as follows: Within the algorithm 
based on Feynman parametrisation the definition of the integration contours is straightforward, 
even in the case of massive propagators. After Feynman parametrisation the singularities in loop 
momentum space are for fixed external momenta and fixed Feynman parameters on a single 
cone. The deformation for the Feynman parameters follows from a kinematical matrix, which 
depends only on the external invariants. However, all terms including integrable singularities are 
raised to the power n, being the number of external particles. This is an artefact of Feynman 
integration. For non-degenerate external kinematics we can have in four space-time dimensions 
maximally four propagators which go simultaneously on-shell. The additional integration over 
the Feynman parameters will effectively lower the degree of the singularities down to the correct 
one. However, this integration is done numerically and is therefore a source of large Monte Carlo 
errors. 

Within the algorithm, which works directly in loop momentum space, the deformation for the 
integration contour is more involved and currently only known in the massless case, but it has the 
advantage that this algorithm does not artificially enhance the degree of a singularity. Within this 
algorithm the singularities lie on the cones (k — qi) 2 = with origins given by qo, q\, q n -\. 
Since 

Pi = qi-qi-u (36) 

the external momenta p, connect the origins of the cones and we arrive for a generic primitive 
one-loop amplitude at the graphical representation shown in fig.® in diagram (a). This is the 
graphical representation of a Wilson loop. For two initial state particles we always have two 
strands in the positive /-direction. This is shown in diagram (a) where we have one strand from 
qo to qj-\ and another strand from qj to q n -i- In the case, where the two initial state momenta 
are adjacent, the diagram degenerates to the one shown in diagram (b) in fig. (QQ). If in addition, 
particle (n — 1) and particle n always couple through an intermediate particle (e.g. photon or 
Z-boson) to the loop, then the propagator corresponding to q n -\ is absent in the loop and we 
obtain the situation shown in diagram (c) in fig. (OQ) consisting of a single strand. 
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t t t 




Figure 1 : A sketch of the loop momentum space: Diagram (a) corresponds to a generic primitive 
amplitude, diagram [b) corresponds to a primitive amplitude, where the two incoming particles 
are adjacent, diagram (c) corresponds to the situation in electron-positron annihilation, where 
the poles due to q n -\ are absent. 



In loop momentum space the points qo, q\, q n -\ are contained within a finite region, which 
we call the "interior region". In diagram (c) of fig. [T] we may define the interior region by the 
intersection of the interior of the backward light cone from q n -i with the interior of the forward 
light cone from qo. The complement to the interior region we call the "exterior region". 

Both algorithms yield identical results, with the direct loop momentum deformation method 
being more efficient for processes with a large number of external particles. The algorithm based 
on Feynman parameters has already been described together with optimisation techniques in 
ref. B2). In this paper we focus on the direct loop momentum deformation method and discuss in 
the next section optimisation techniques for this method. 

3 Optimisation techniques 

In this section we discuss the employed optimisation techniques for the Monte Carlo integration 
over the loop momentum space and the phase space of the final state particles. The contour of 
integration is deformed directly in loop momentum space, without the introduction of Feynman 
parameters. The optimisation techniques are a combination of standard Monte Carlo optimisation 
techniques H2l (which can be applied to any Monte Carlo integration and do not require any 
particular information on the integrand), and improvements obtained by taking into account the 
physical nature of the problem. We point out that in the latter case we only use information which 
follows from the general principles of quantum field theory, for example that the singularities of 
the integrand are given by the poles of the propagators. We do not use any process-dependent 
information. Therefore, all optimisation techniques are process-independent. 
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3.1 Holomorphic and non-holomorphic division into sub-channels 

We start with a general observation. Consider a complex contour integral 

/ = Jdzf(z), (37) 
c 

where f(z) is a meromorphic function. Suppose f(z) = fi(z) + fi(z), where f\{z) and f 2 (z) are 
again meromorphic functions of z. Then 

/ = J dzMz) + J dzf 2 (z). (38) 
Ci c 2 

The contours C\ and 62 may differ from the original contour C, as long as no poles are crossed 
in going from C to C\ or C 2 - We can use this fact by optimising the contours for f\ and f% 
separately. 

Suppose on the other hand that f(z) = gi(z,Z*) + g 2 (z,Z*), where gi(z,z*) and g 2 (z,z*) are 
now - taken individually - non-meromorphic functions. If we now divide the integral into two 
channels, 

/ = Jdzgi(z,z*)+fdzg 2 (z,z*), (39) 
c c 

we have to use the same contour for both channels. However, we may use different parametri- 
sations of the same contour. This is the typical situation when we set gi(z,z*) = wi(z,z*)f(z), 
g2(z 7 z*) = W2(z,z*)f(z) with weight functions w\{z,z*) and W2(z,z*), satisfying w\ +w 2 = 1 
and involving the complex modulus. This will be discussed in more detail in section l3".4.2l 
Let us now return to the integral in eq. (|32l : 



/ 




We define 

/uv« = Ud—^- (41) 

j=0 K 

Clearly, fuv(k) is a meromorphic function of k, with poles only at k — fijjy = 0. We can use the 
function fuv(k) to split our original integral / into two channels 

/ = hxt+hnu (42) 
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with 



'ext 



'int 



d 4 k 



d 4 k 



fw(k) 



R(k) 



n—2 / 

j=o v 



2 
j 



[Wuv«] 



R{k) 



n-2 

n 

7=0 



m 



(43) 



This splitting is holomorphic in k, therefore we can evaluate / ex t and I mt with two different con- 
tours. The integrand of 7 ext has only poles at k — A'uv = ant ^ can ^ e evaluated with a relatively 
simple contour. Since k = k — Q we further have within I mt 



WuvW 



1 



n-2. 



,n-l 



If we define the arbitrary four- vector Q to be 

Q 



(k 2 )" 2k 




1)2 



+ 0(k 



2^^-4^ 



.(44) 



n-2 



(45) 



we can arrange that the integrand of 7i nt drops off with two extra powers of k for \k\ — > °°. This 
ensures that /j nt receives less contributions from the ultraviolet region. 

In summary, the division into the two channels of eq. (1421) has the advantage that the integral 
Text has a simple pole structure, while the integral 7i nt drops off with two additional powers in the 
ultraviolet region. 



3.2 Improvement of the ultraviolet subtraction terms 

Ultraviolet finiteness requires that the integrand falls of stronger than 1 /\k\ 4 for \k\ — > °°. The 
ultraviolet subtraction terms from ref. [0 ensure that the integrand falls of like 1 /\k\ 5 for \k\ — > 
oo. This ensures that the integral is ultraviolet finite. However, it turns out that the ultraviolet 
(or exterior) region contributes significantly to the statistical error. In an expansion around the 
ultraviolet propagator (k — A^uv) 1 a term which falls off like 1 /\k\ 5 for \k\ — > °° is given by 

d>k hX (46) 



where X is a four-vector independent of k. This term integrates to zero. At the integrand level 
this term is oscillating. The integrand changes the sign under k — > (—k). Oscillating terms pose 
no problem for a Monte Carlo integration, if the amplitude of the oscillation is small compared 
to the terms which give a non-vanishing contribution after integration. Unfortunately it turns out 
that the amplitude of the oscillations related to ultraviolet terms of order 1 /\k\ 5 or 1 /\k\ 6 is not 
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2 3 4 5 6 




Figure 2: The ultraviolet behaviour of some example diagrams. To the right the number of 
external particles increases, in the vertical we have the various powers of the large \k\ -behaviour, 
ranging from \k\~ 2 (top) to \k\~ 6 (bottom). 



small. The amplitude is enhanced whenever an external invariant approaches the jet resolution 
parameter. 

The situation can be improved by a better damping in the ultraviolet region. To this aim we 
modify the ultraviolet subtraction terms such that they subtract out also the l/|fc| 5 - and the l/|fc| 6 - 
behaviour of the integrand, whenever there is a possibility of having a small two-particle external 
invariant. This amounts to modifying the ultraviolet subtraction terms for the propagators and 
the three-valent vertices, as well as the soft and collinear subtraction terms. 

Note that if one would like to achieve that the complete integrand has a behaviour better than 
l/\k\ 6 in the ultraviolet region, then it is not sufficient to include to the modifications mentioned 
above the corresponding one for the four-gluon vertex. In addition one has to include at 0(|fc|~ 6 ) 
subtraction terms for new four-, five- and six-valent vertices. This is illustrated in fig. (0, where 
we show the ultraviolet behaviour for some example diagrams. From a technical point of view 
the inclusion of new subtraction terms for four-, five- and six-valent vertices is not a problem, but 
the inclusion would modify the recursion relations for the ultraviolet subtraction terms. There is 
a trade off between the level of improvement in the subtraction terms and the efficiency for the 
numerical evaluation of the integrand. We have chosen to improve only the two- and three-valent 
parts of the ultraviolet behaviour. Empirically it turns out that this is sufficient to reduce the 
oscillating behaviour in the ultraviolet region. 
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Technically the improvement is done as follows: We have to improve eq. (|29l to order \k\ 1 . 
This follows easily from 



1 



1 



(k — p) 2 — m 2 
and the expansion 



k 2 - 2 



A'uv 



1 - 



2k-(p-Ql + (p-Qy-m 2 +i J 2 ]W 



k 2 - 



WJV 



k 2 



(47) 



[l-ax + bx 2 ] 1 = l+ax+ (a 2 -b)x 2 + (a 3 -2ab)x 3 + (a 4 -3a 2 b + b 2 ) x 4 



with 



2k-(p-Q) 
k 2 



(p-Q) 2 -m 2 + iu 2 JY 



k 2 - 



(48) 



(49) 



We apply this expansion to the propagator and the three- valent vertex corrections and determine 
the ultraviolet behaviour up to order \k\~ 6 . We then add appropriate terms of order |fc|~ 8 (which 
are beyond the order to which we are working) to ensure that the integral over the ultraviolet 
subtraction terms is not changed. In other words, we modify only the unintegrated version of the 
ultraviolet subtraction terms in such a way that the integrated version of the ultraviolet subtrac- 
tion terms remains unchanged. Leaving the integrated version of the ultraviolet subtraction terms 
unchanged ensures that the sum of all ultraviolet subtraction terms is proportional to a Born am- 
plitude. The explicit forms of the improved ultraviolet subtraction terms are rather lengthy. We 
therefore do not list them here, however they can be obtained systematically with the procedure 
outlined above. 

In addition we improve the ultraviolet behaviour of the soft and collinear subtraction terms, 
such that they fall off like 0(|fc|~ 7 ) in the UV-region. We give here the improved subtraction 
terms for massless QCD, the extension to massive quarks is straightforward. The improved soft 
subtraction term is given by 



soft 



1 E *Pj 



Pj+l 



1 



1 



,(°) 



(50) 



The second term in the square bracket is new and subtracts off the leading |fc|~ 6 -behaviour. 
Note that an individual soft subtraction term is proportional to a Born amplitude. Therefore we 
can easily allow for a modification of the integrated soft subtraction term. Integrating the soft 
subtraction term in eq. (|50l we obtain 



V 8 



/ 



d u k (1) 
{2n) D soft 



1 



(47t) 2 r(i 

+ 0(8). 



~ 2 PJ-Pj+1 

P 2 



+ 



2 PJ-Pj+1 

A'UV 



(51) 
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The improved collinear subtraction is given by 



.(1) 

'coll 



*E(-2) 



S jgvw (k)_ v kfj S j+ 1 guv (k 2 j ,kj +l ^ 



1(0) 



k 2 k 2 



k 2 k 2 



(52) 



where the function guv reads 
guv (kj_ h kj) = 1 



k 2 k 2 



k 2 _ i k 2 (2k-q j _ l + 2k-q j ) 



[{k-QY-^f 



[(k-Q) 2 -vhv} 3 
(2k ■ qj- i) 2 +{2k- qj) z + (2k- q^) (2k ■ qj) 



k 2 k 2 



(53) 



- \2 



[(k-Q) 2 -^) 3 
In eq. (1531) we have used the notation 

9j- i = 4]-i-Q, q~j = <ij-Q- 

Integrating the collinear subtraction term we obtain 



[(k-Q) 2 -^]' 



(54) 



-1..2e 



d D k m 



(2%) 



D coll 



(4jc)2r(l-e)^ 



Af + 0(e). 



(55) 



Let us summarise what we achieved so far: We have split the integral in eq. (1401 ) into two sub- 
integrals 7 ext and I mt . With the improved ultraviolet subtraction terms, the integrand of 7 ext falls 
off like \k\ ~ 7 for all parts corresponding to one-loop rc-point functions for n < 3. The integrand 
of hxt fells off like \k\~ 5 for the parts corresponding to one-loop n -point functions for n > 4. On 
the other hand, the integrand of 7i nt is suppressed by two additional powers of \k\. Therefore the 
parts corresponding to one-loop n-point functions with n < 3 in 7j nt fall off like \k\~ 9 , while the 
parts corresponding to one-loop n-point functions with n > 4 fall off like \k\~ 7 . 

We can further improve the ultraviolet behaviour of all terms by another power of \k\ as 
follows: All terms which scale with an odd power in the ultraviolet region are necessarily anti- 
symmetric under the substitution k — > (—k). We can eliminate these terms by sampling simulta- 
neously the points k and (—k). With this method we can reduce the leading ultraviolet behaviour 
(which scales like an odd power in all cases above) to the next lower even power. This will be 
discussed in more detail in sub-section l3.4l 



3.3 Optimisation of the integration contour 

Since the division of the integral in eq. (|40l into the two sub-integrals 7 ext an d 7i nt is holomorphic, 
we may use different integration contours for 7 ex t and 7^. 
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3.3.1 The contour deformation for Text 



By construction the integrand of 4xt has only poles at k 2 — = 0. The poles are located on a 
single light cone in loop momentum space and the contour deformation is rather simple. For 7 ext 
we use the contour deformation 

k = k + iK (56) 



with 



(57) 



(The position of the indices is correct. We use the metric tensor g^ v = diag( 1,-1,-1,-1) to 
indicate that the spatial components of contravariant four- vector are the negative of the spatial 
components of the contravariant four- vector W 1 — Q^.) We have 



72 2 

k -a/uv 



2i(k-Q)o(k-Q) 



2 

A'UV' 



(58) 



where aob denotes the Euclidean scalar product of the four- vectors a and b. As already men- 
tioned we will always choose /j^y purely imaginary with Im /j^y < 0. Therefore the term — ^uv 
alone ensures that the imaginary part of k 2 — /j^y is always positive. The reader may ask why 
we deform the integration contour at all. The answer is given by the ultraviolet behaviour of the 
propagator [k 2 — fi^y]^ 1 : With the deformation as in eq. (|56T ) the propagator falls off always like 
\k\ " 2 for \k\ — > oo, while with no deformation the propagator would remain constant along the 
light cone k = 0. The Jacobian of the contour deformation is given by 



dkf 



-4i. 



(59) 



3.3.2 The contour deformation for /j nt 

The integrand of Ii ni has good ultraviolet properties, but it has a more complicated infrared struc- 
ture. With fj^jy purely imaginary and Im /j^y < 0, we have to take into account the poles given 
by 

k)-m 2 = 0, 0<j'<n-2. (60) 

In this section we discuss the massless case mj — 0, so the poles are given by k 2 = for < j < 
(n — 2). For I m we use the contour deformation along the lines of ref. [6]. We first define 

P = ^(qo + qn-i). (61) 

The four-vector P is the centre of the forward light cone from qQ intersected with the backward 
light cone from q n -2- We recall that the four-vector Q has been defined in eq. (l45l) as the average 
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of the (n— 1) four-vectors qo, q n -2- Note that in general P and Q are not equal. We further 
define 



2fc;;-2 • Pn- l 

(q n -2-qo) 



2 ' 



2^/7— 1 ■ Pn 
(qn-2-qo) 



2 ■ 



We set 



n-2 



k = k + i\x, = (c+ - c-) {q n -2 - qo) - £ c jh- 



The coefficients c± and c/ (0 < j < n — 2) are given by 



c+ = (x + x)e(x+x)h-(k n - 2 )g-(k-P), 

c- = (— x— x) 9 (— x — x)h+ (ko) g+ (k — P) , 

c 7 = h + (k j+l )h^(kj^)g(k-P), l<j<n-3, 

: h-(k n -3)g{k-P). 



C n -2 



The functions /z + , g + , g_ and g are defined in [6] and read 



(62) 



(63) 



(64) 



h+ (k) 



k°) +M] 



h-{k) 



+ k°) +M] 



+ k° ), (65) 



J) M \,v 8± (k) = (66) 



kok + MV 



1+|1± *° 



These functions depend on the five parameters M\, M%, M3, yi and y 2 . The default choice for 
these parameters is 



M 1= 0.05J^(^-2-^o) 2 , M 2 =M 3 = J^(q n -2-qo) 2 , (67) 



yi=0.7, y 2 = l. (68) 

These default parameters have been given in ref . [6] . We have varied these parameters in different 
processes. It turns out that these values are a good choice for a wide range of processes. It 
remains to define the scaling parameter X. We set 

1 = min(l,X ,...A/7-2AuvAcoii), (69) 
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with Xj (0 < j < n — 2) given by 



K 2 k 2 j 



for 0<2(k-^) 2 <k 2 ^, 
for 0< k 2 £2 < 2(k-£ 7 -) 2 , 
for K 2 fc 2 <0<2(K-y 2 . 



4(T^-igg 

(2K 2 ) 2 
4(K-^) 2 -2K 2 fc 



(70) 



A-coli is given by 



Finally, Xuv is given by 



coll 



1 



1 

4C' 



n-2 
7=0 



J- 



(71) 



4(*-e)-K 



for 4 (jc — Q^j ■ K > Inijtf 
otherwise. 



uv 



(72) 



Xq, Ki-2 and X-con are defined as in ref. [|6)- In addition we introduced which protects 
the ultraviolet propagator lc — fi^y from going on-shell due to a too large contour deformation. 
Note that X\jy differs from 1 only for 



4(^-2)-K<Im A /f JV <0. 



(73) 



It can be shown that with the contour deformation as in eq. (1631) the propagators fall off like 
\k\ ~ 2 for \k\ — > °°. For the proof it is convenient to discuss the time-like, space-like and light-like 

regions separately. In the time-like (k — > °°, k = const) and space-like (\k\ — > °°, k = const) 
regions the contour deformation goes to zero in the limit of \k\ — > °°. Therefore k 2 « k 2 and 
the large |£|-behaviour is given by the real part. In the light-like region (Jc 2 ~ 0) the contour 
deformation does not vanish in the limit where \k\ — >■ °° and the imaginary part of k 2 scales like 
\k\ 2 . 

The Jacobian 



dkf 



dk v 

for the contour deformation for 7i nt is computed numerically. 

3.4 Sampling in loop momentum space 

After the contour deformation we have a four-dimensional integral 

d 4 k 



(74) 



(2%) 



lf(k) 



(75) 
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over the four real variables k°, k l , k 2 and P. For each dimension the integration is from minus 
infinity to plus infinity. In order to be able to use standard methods for the Monte Carlo inte- 
gration like Vegas [|43ll44| we map M 4 to the four-dimensional unit hyper-cube [0, l] 4 . There 
are many possible choices for such a mapping. We choose a mapping which approximates the 
peak structure of the integrand. This amounts to importance sampling and reduces the statistical 
Monte Carlo error. In order to further minimise the Monte Carlo integration error we use where 
possible the method of antithetic variates. Within this method we always evaluate a few points 
together, which are anti-correlated. In this way oscillations are significantly reduced. 



3.4.1 Sampling for I t 



cxl 



We start the discussion with our method for sampling the integrand of Text- We generate the real 
vector k iea \ =k — Qas follows: We use four uniformly in [0, 1] distributed random numbers uq, 
and define the four quantities kg, and § by the equations 



tan | M , 

£, = arccos (1 — 2«i) 
= arccos (1 — 2«2) 
§ = 2nu3. 



(76) 



/ji is an arbitrary scale, which we take to be of the order of the centre-of-mass energy. We then 
set 



sin sin (|) 

^real = focOS^, K = ^-E sin£, k Ka \ = k r [ sin0COS(|) 

COS0 



(77) 



The Jacobian of this transformation is 



dk 



du 



2*4(4+A<?)si< 



(78) 



In section [3T2l we have shown that the integrand of 7 ext falls off like |/:|~ 7 for all parts corre- 
sponding to one-loop n-point functions for n < 3, and like |&|~ 5 for the parts corresponding to 
one-loop n-point functions for n > 4. By sampling always the two points with loop momenta k rea \ 
and (— ^reai) together we can reduce the ultraviolet behaviour to \k\~ s and \k\~ 6 , respectively. 



3.4.2 Division into sub-channels for 7i nt 

The efficient sampling of the integrand of 7i nt is more involved. We recall from fig. (Q3 that the 
integrand is characterised by a strand of (n — 2) line segments (diagram (c) of fig. ©). We 
view a single line segment as a basic building block and we divide 7i nt into (n — 2) sub-channels, 
such that each sub-channel corresponds to a line segment. This is shown pictorially in fig. ©. 
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Figure 3: Diagram (a) shows the origins of the light cones in electron-positron annihilation. The 
origins are given by the (n— 1) vertices. The vertices are connected by (n — 2) line segments. 
We decompose 7i nt into (n — 2) sub-channels, such that each sub-channel corresponds to one line 
segment. This is shown in diagram (b), where the line segment from q\ to q% is drawn. 



Technically this is done as follows: After contour deformation we have an integral over a real 
four-dimensional space 



'int 



d 4 k 

j2%r 



The line segments correspond to 

k = qj+x(qj+i—qj) , 0<x<l, 0<j<n — 3 



(79) 



(80) 



in loop momentum space. We have qj + \ —qj = Pj+\- We split the original integral into several 
channels, such that each critical line segment corresponds to a separate channel. We then use for 
each channel a dedicated mapping for the loop momentum. The splitting into different channels 
is done as follows: We re- write the original integral as 



'int 



n-3 

E 

!=0 



(*)/(*) 



with 



n-3 



W 



i>0 and £w;(&) = l. 

(=0 



(81) 



(82) 



Since the sum over all weights W/ equals one, the sum over all channels equals the original 
integral. For the weights wi we use 

a 



W, 




a > 



(83) 



k 2 \\k 2 
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with a = 2. Note that for the weights we take the norm of the complex quantities kj. Therefore 
the division into the sub-channels is not holomorphic and we have to use the same deformation 
contour for all sub-channels. The weights have the properties that 



lim wi = 1 if k ->■ qi + x(q i+ i- qi 



and 



lim Wi = if £— >• qj+x(qj + i—qj) with i ^ j. 
This ensures that in each channel there is only one critical line segment. 



(84) 



(85) 



3.4.3 Sampling of an individual sub-channel of I\ 



int 



After having divided I mt into different sub-channels we can discuss how to sample the integration 
points for a specific sub-channel. Recall that a line segment in loop momentum space goes from 
qj to qj + \ . The four- vector from qj to qj + i is given by pj + \ =qj + \ — qj. Since we are discussing 
a single line segment we will in the following simply write p instead of pj+i- For massless 
external particles p is a light-like four- vector. We now look for an appropriate coordinate system 
in the case when there is a distinguished vector p. A possible choice are generalisations of 
elliptical or prolate spheroidal coordinate systems to four dimensions. 

In detail we generate k as follows: Let p be a light-like four- vector. In spherical coordinates 
p can be written as 

p° = |p|cos0i, 

p l = \p\ sin0icos02, 

p 2 = \p\ sin0i sin02cos(|)3, 

p 3 = \p\ sin0i sin02sin(|)3, (86) 



where 



arccos 



02 = arctan 



(P 2 ) +(P 2 



§3 



arctan 



(87) 



The angles 0i, 02 and §3 define three rotation matrices 



#1 



/ COS 01 

sin 0i 


\ 



— sin 0i 
cos 01 







*3 






\ 




( 1 ° 










COS 02 


1 





sin 2 





1 ) 




\0 


f 1 








\ 







1 
















COS(])3 — 


sin(j)3 









sin(])3 cos(|)3 J 





\ 

-sin 2 

COS 02 

1 J 



(88) 



(89) 
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We generate the (real) loop momentum k as follows: 



where p = pj + \ and 



q j + ^p + R 3 -R 2 -Ri-~k' 1 



( coshp cos£, 



(90) 



k' 



\ 



sinhp sin£, cosG 
sinhp sin£, sinG cos(|) 
y sinhp sin£, sinG sin(]> J 



(91) 



p, ^, G and § are generalised elliptical coordinates. The coordinate p ranges from [0,°°[, the 
coordinates % and G have the range [0,7t], whereas § is in the range [0, 2ft], The Jacobian is given 
by 



dk 



dk' 



|/?| 4 sinh 2 p sin 2 £, sinG (sinh 2 p + sin 2 ^) . 



(92) 



The variables (p, ^, 0, <f») we generate as follows: 

, A uq 71 

p = In l + i— r tan -mo 
V P 2 



G 




(l+B)(^)- 2 " 2 -8 

e-(l+e)(i?) 



8 



< u 2 < j, 

j <U2 < 1, 



with 



sinhp sin^. 



(93) 



(94) 



jUo is again an arbitrary scale, which we also take to be of the order of the centre-of-mass energy. 
The functions in eq. (1931) have been chosen such that they approximate the typical peak structure 
of the integrand. The Jacobians are 



duo 
du i 

ae 

d"3 



2 
71, 



4 + (eP-l) 2 

\p\ 



(£+|cosG|) 



sinG 



In 



1+e 



2%. 



(95) 
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Figure 4: The recurrence relation for the pure gluon current . 



Again we use the method of antithetic variates: We observe that the integrand has a periodic 
behaviour in (|). Therefore combining the evaluations at (|) and ((|) + 7t) mod (2n) averages out these 
oscillations. In addition, the integrand is for p — > strongly peaked and antisymmetric around 
= 7t/2. Evaluating the integrand at and 71 — averages out this behaviour. Furthermore 
we evaluate the integrand always with the values k' and (—k) in eq. ( |90l . This improves the 
ultraviolet behaviour. 

4 Recurrence relations 

This section is devoted to the computation of the integrands. The computation of the integrands 
is done efficiently with the help of recurrence relations. We first discuss in sub-section 14.1 1 tree- 
level recurrence relations, which are directly relevant to the Born contribution, the integrated 
subtraction terms and the unintegrated soft and collinear subtraction terms. In sub-section 14.21 
we discuss how the integrand of the bare one-loop amplitude is calculated. Cutting the loop open 
reduces the computation to a tree-like problem. In sub-section !4.3l we treat the computation of the 
ultraviolet subtraction term. Again, this is done recursively with a modified tree-like recurrence 
relation. 

4.1 Tree level recurrence relations 

We first review the computation of tree-level partial amplitudes. Berends-Giele type recurrence 
relations [431 build tree-level partial amplitudes from smaller building blocks, called colour or- 
dered off-shell currents. Off-shell currents are objects with n on-shell legs and one additional 
leg off-shell. Momentum conservation is satisfied. It should be noted that off-shell currents are 
not gauge invariant objects. Recurrence relations relate off-shell currents with n legs to off-shell 
currents with fewer legs. The recursion relation for the pure gluon off-shell current is depicted 
in figlH The recursion starts with the one-currents: 

jf\l) = $(pi,qi)- (96) 

is the polarisation vector of the gluon corresponding to the polarisation X, pi the four-momen- 
tum of the gluon and qj an arbitrary light-like reference momentum used to define the polarisation 
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of the gluon. The recursive relation states that in the pure gluon off-shell current a gluon couples 
to other gluons only via the three- or four-gluon vertices: 



j£\m,...,n) 



where 



p2 



n-l 



£<^-P M) n,P wJ ,P; + l,„)/f (m,...,j)7^(7 + l,...,n) 



r(«), 



j=m 



n—2 n—1 



+ [ I Vr C 4°\m,...J)4°\j+l,...M 0) (l + h..,n) 

j=m l=j+l 



(97) 



Pij = Pi + Pi+i+...+Pj 
and V3 and V4 are the colour ordered three-gluon and four-gluon vertices 



(98) 



V, 



/jvpa 



(99) 



From an off- shell current one easily recovers the on- shell amplitude by removing the extra prop- 
agator, taking the leg (n+1) on-shell and contracting with the appropriate polarisation vector: 



A(°)(l,..,n+1) = e^ Pn+h q n+1 )iPl n 4 0) (l,...,n„ 



l,n = ~ Pn+\ 



(100) 



Similar recurrence relations can be written down for the quark and antiquark currents, as well 
as the gluon currents in full QCD. If there is only one quark line, we have for the quark and 
antiquark off- shell currents 



C/ (0) (m,...,n) 



V {0) (m,...,n) 



n-l 



m,n 



i£t/ w (m,...,0V^(/ + l,...,n)-^ 

i=m r m,,i 

-i^E^(«,-",OV (0) (/+l,...,»). (101) 



P 2 

1 m,n i= m 



This is shown pictorially in figJ51 The recursions start with U^°\l) = u(pi) and V^°\l) = v(pi), 
respectively. The colour ordered quark-gluon vertex is given by 



V* 1 - 

qgq 



-if. 



(102) 



Note that in eq. 

(USD the quantities U {0) , V^, V qg9 and If are matrices in Dirac space, and the 
order is relevant. The partial amplitudes with one quark-antiquark line are given by 



A(°)(l 9 ,2,...,n-1,^ 



™(piW 2n V(°)(2,..,n-l,n) 



(103) 
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Figure 5: The recurrence relations for the quark current and the antiquark current. 

The quark and antiquark off- shell currents enter also the computation for the amplitudes 

e + e~ ->■ q,g,..., g,q. (104) 
Taking all particles as outgoing we have 

A(°)(l„ 2 g , (n - 3) g , (n - 2) 9 , (n - \) h nj) = 

£Z7 (0) (i,...,oy^ Y/z y (0 H^i,--^-2)^ E >-i,n). (io5) 
(=i 

Here, V^- ^ z is the electroweak quark-photon/Z-boson vertex and J^ w (n — 1 ,n) the electroweak 
current, including the photon and the Z-boson propagator. 

4.2 One-loop recurrence relations 

We now turn to the computation of the integrand of the bare one-loop amplitude. For a given loop 
momentum k and external momenta p\, p n we consider the unintegrated one-loop currents 

y^ 1 ) (unintegrated one-loop pure gluon current), (unintegrated one-loop quark current) and 
V^ 1 ) (unintegrated one-loop antiquark current). Note that for a given k and given p\, p n 
all momenta are fixed, in particular kj = k — qj, qj = p\ + ... + p,-. The recursion relations 
for the unintegrated one-loop currents [14611471 are shown in fig. [6] for the gluon current and in 
fig- 13 for the quark and the antiquark current. In the recurrence relations for the unintegrated 
one-loop currents there are two types of vertices. For the first type the off- shell leg couples to 
a tree-like vertex, in which case the recursion relates the unintegrated one-loop off-shell current 
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Figure 6: The recurrence relation for the unintegrated one- loop gluon current. 



n+ 1 




n+ 1 




Figure 7: The recurrence relation for the unintegrated one- loop quark current and the uninte- 
grated one-loop antiquark current. 
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Figure 8: By cutting open the loop at a gluon line one replaces the tensor structure of the indicated 
gluon loop propagator by a sum over four (pseudo-) polarisations. 



to an unintegrated one-loop off-shell current with fewer legs. This recursion terminates with 
unintegrated one-loop one-currents, which are zero. This corresponds to the fact that self-energy 
corrections on external lines are not included in the computation of the integrand of the bare 
one-loop amplitude. 

For the second type the off- shell leg couples directly through a vertex to the loop. We call 
this contribution the "direct contribution". In this case two edges of the vertex are connected to 
loop propagators. We can cut open one of these two edges by replacing the tensor structure of 
the corresponding propagator by a sum over (pseudo-) polarisations. Technically, this is done 
as follows: If the edge we would like to cut open corresponds to a gluon, we replace in the 
numerator of the gluon propagator in Feynman gauge by 

4 

Sa* = E4 4 s (106) 
(=1 

with four standard (pseudo-) polarisations 

4 1} = (1,0,0,0), 4 2) = (0,-/,0,0), 
4 3) = (0,0,-i,0), 4 4) = (0,0,0,-/). (107) 

This is shown in fig. [81 If, on the other hand, the edge corresponds to a (massless) quark line, we 
replace in the numerator of the quark propagator by 

H = t + ^-^ k 9 = k-^- qi (108) 
2kq 2kq 

where q is a light-like reference momentum and k 9 is by construction light-like. We then replace 
ly and q 1 by a polarisation sum 

f = £ u(k\X)u(k\X), 

4 = Y,u{q,X)u(q,X). (109) 

This is shown in fig.|9] The construction above for massless quarks can be generalised to massive 
quarks. 
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Figure 9: By cutting open the loop at a quark line one replaces the tensor structure of the indicated 
quark loop propagator by a sum over four (pseudo-) polarisations. 




Figure 10: The recurrence relation for the the pure gluon chain. The chain results from cutting 
open the loop. 

In the case of a one-loop gluon current we also have to cut open the ghost loop. Since the 
ghosts are scalar particles, this is rather simple and does not involve a non-trivial polarisation 
sum. 

As a sideremark we note that within the context of calculations based on Feynman diagrams 
the technique of cutting open the loop has recently been discussed in ref. Il48l . 

In all cases after cutting open the loop we obtain an object which we call a chain. This object 
can again be calculated recursively. The recursion relation corresponds to a tree-like calculation 
with (n+l) external legs. One of the external legs corresponds to the cut loop propagator. 
This leg is not on-shell. However this does not affect the recurrence relations. We show 
two examples for the recurrence relations of chains. In fig. [T0]the recurrence relation for the 




Figure 11: The recurrence relation for the the clockwise ghost- antighost chain. The chain results 
from cutting open the ghost loop. 
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Figure 12: The recurrence relation for the unintegrated one-loop electroweak quark-gluon- 
antiquark current. 

pure gluon chain is given. As a second example we show in fig. [TT]the recurrence relation for 
the clockwise ghost- antighost chain. In addition there are further chains: From cutting open all 
direct contributions to the unintegrated one-loop gluon current, we also obtain the anti-clockwise 
ghost- antighost chain. Cutting the direct contribution to the unintegrated one-loop quark current 
gives a quark-gluon chain, cutting the corresponding contribution to the unintegrated one-loop 
antiquark current gives a gluon-antiquark chain. These chains have similar recurrence relations, 
which we do not show explicitly. 

For completeness we mention that for the process e + e~ —¥ q,g,...,g,q we also need the un- 
integrated one-loop electroweak quark-gluon-antiquark current. This current can be calculated 
with the same methods as discussed above. The recurrence relation for this current is shown in 
fig- ELU In the direct contribution we cut open the quark propagator as outlined in eq. (11081 ) and 
eq. (1 1 09b . We obtain an antiquark-quark chain, which again can be computed with a tree-like 
recurrence relation. 

4.3 UV recurrence relations 

Let us now turn to the computation of the ultraviolet subtraction term. We recall that the basic 
building blocks of the ultraviolet subtraction term are propagator and vertex subtraction terms. 
We can think of the ultraviolet subtraction term as a sum over diagrams, where each diagram 
has a tree-structure with exactly one propagator or vertex replaced by a basic ultraviolet subtrac- 
tion term. The complete ultraviolet subtraction term can again be calculated recursively. The 
recursion relation is shown in fig. [13] for the gluon current and in fig. [14] for the quark current. 
The recurrence relation for the antiquark current has a similar structure as the recurrence rela- 
tion for the quark current and is not shown explicitly. The structure of the recursion relations 
is as follows: Either the off-shell leg is attached through a Born propagator to a Born vertex. 
In this case exactly one of the sub-currents carries a basic ultraviolet subtraction term, while all 
other sub-currents are Born currents. Or the off-shell leg is connected through a basic ultravio- 
let subtraction term to the sub-currents. The basic ultraviolet subtractions terms can be either a 
propagator subtraction term or a vertex subtraction term. In this case all sub-currents are Born 
currents. 

For the process e + e~ — > q, g, . . . , g, q we also need the off-shell current including the ultraviolet 
subtraction term for the electroweak quark-gluon-antiquark vertex. The corresponding recursion 
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Figure 13: The recurrence relation for the ultraviolet subtraction term to the pure gluon current. 



n+1 





Figure 14: The recurrence relation for the ultraviolet subtraction term to the quark current. 
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Figure 15: The recurrence relation for the ultraviolet subtraction term to the electroweak quark- 
gluon-antiquark current. 

relation is shown in fig. [151 The basic ultraviolet subtraction term for the electroweak quark- 
gluon-antiquark vertex can be obtained from the sub-leading N c quark-gluon vertex subtraction 
term by adjusting the coupling factor appropriately. 

5 Conclusions 

In this paper we have given a detailed account on the techniques used to improve the efficiency 
of the Monte Carlo integration in a numerical approach for the computation of one-loop QCD 
amplitudes. The techniques fall into two categories: Techniques in the first category reduce the 
statistical error of the Monte Carlo integration. This is done by dividing the integration into 
sub-channels and optimising the integration in each sub-channel. Of particular importance is the 
improvement of the ultraviolet subtraction beyond the formally required \k\~ 5 fall-off. Within 
the second category we discussed CPU-efficient methods for the computation of the integrands. 
This is done with the help of recurrence relations. 

All techniques discussed in this paper are process-independent or have an obvious generali- 
sation towards more complex processes. 
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